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ABSTRACT 

Multiplicative noise models occur in the study of several coherent 
imaging systems, such as synthetic aperture radar and sonar, and ul- 
trasound and laser imaging. This type of noise is also commonly 
referred to as speckle. Multiplicative noise introduces two additional 
layers of difficulties with respect to the popular Gaussian additive 
noise model: (1) the noise is multiplied by (rather than added to) 
the original image, and (2) the noise is not Gaussian, with Rayleigh 
and Gamma being commonly used densities. These two features 
of the multiplicative noise model preclude the direct application of 
state-of-the-art restoration methods, such as those based on the com- 
bination of total variation or wavelet-based regularization with a 
quadratic observation term. In this paper, we tackle these difficul- 
ties by: (1) using the common trick of converting the multiplicative 
model into an additive one by taking logarithms, and (2) adopting the 
recently proposed split Bregman approach to estimate the underlying 
image under total variation regularization. This approach is based on 
formulating a constrained problem equivalent to the original uncon- 
strained one, which is then solved using Bregman iterations (equiv- 
alently, an augmented Lagrangian method). A set of experiments 
show that the proposed method yields state-of-the-art results. 

Index Terms — Speckle, multiplicative noise, total variation, 
Bregman iterations, augmented Lagrangian, synthetic aperture radar. 

1. INTRODUCTION 
1.1. Coherent Imaging and Speckle Noise 

The standard statistical models of coherent imaging systems, such as 
synthetic aperture radar/sonar (SAR/SAS), ultrasound imaging, and 
laser imaging, are supported on multiplicative noise mechanisms. 
With respect to a given resolution cell of the imaging device, a co- 
herent system acquires the so-called in-phase and quadrature com- 
ponent^ which are collected in a complex reflectivity (with the in- 
phase and quadrature components corresponding to the real and imag- 
inary parts, respectively). The complex reflectivity of a given resolu- 
tion cell results from the contributions of all the individual scatterers 
present in that cell, which interfere in a destructive or constructive 
manner, according to their spatial configuration. When this configu- 
ration is random, it yields random fluctuations of the complex reflec- 
tivity, a phenomenon which is termed speckle. The statistical prop- 
erties of speckle have been widely studied and there is a large body 

' The in-phase and quadrature components are the outputs of two demod- 
ulators with respect to, respectively, cos{ujot)) and sin{uJot), where uja is 
the canier angular frequency. 



of literature 1111 . 1141 . Assuming no strong specular reflectors and 
a large number of randomly distributed scatterers in each resolution 
cell (relative to the carrier wavelength), the squared amplitude (in- 
tensity) of the complex reflectivity is exponentially distributed 1 14|. 
The term multiplicative noise is clear from the following observa- 
tion: an exponential random variable can be written as the product 
of its mean value (parameter of interest) by an exponential variable 
of unit mean (noise). The scenario just described, known as fully 
developed speckle, leads to observed intensity images with a charac- 
teristic granular appearance due to the very low signal to noise ratio 
(SNR). Notice that the SNR, defined as the ratio between the squared 
intensity mean and the intensity variance, is equal to one (OdB). 

1.2. Restoration of Specified Images: Previous Work 

A common approach to improving the SNR in coherent imaging con- 
sists in averaging independent observations of the same pixel. In 
SAR/SAS systems, this procedure is called multi-look (A/-look, in 
the case of M looks), and each independent observation may be ob- 
tained by a different segment of the sensor array. For fully developed 
speckle, the SNR of an A/-look image is M. Another way to obtain 
an M-look image is to low pass filter (with a moving average kernel 
with support size M) a 1-look fully developed speckle image, mak- 
ing evident the tradeoff between SNR and spatial resolution. A great 
deal of research has been devoted to developing nonuniform filters 
which average large numbers of pixels in homogeneous regions yet 
avoid smoothing across discontinuities in order to preserve image 
detail/edges 1 9 1 . Many other speckle reduction techniques have been 
proposed; see 1 14| for a comprehensive literature review up to 1998. 

A common assumption is that the underlying reflectivity image 
is piecewise smooth. In image restoration under multiplicative noise, 
this assumption has been formalized using Markov random fields, 
under the Bayesian framework 1141 , |2| and, more recently, using 
total variation (TV) regularization |T|, 1121 . 1161 . |17| . 

1.3. Contribution 

In this paper, we adopt TV regularization. In comparison with the 
canonical additive Gaussian noise model, we face two difficulties: 
the noise is multiplicative; the noise is non-Gaussian, but follows 
Rayleigh or Gamma distributions. We tackle these difficulties by 
first converting the multiplicative model into an additive one (which 
is a common procedure) and then adopting the recently proposed 
split Bregman approach to solve the optimization problem that re- 
sults from adopting a total variation regularization criterion. 

Other works that have very recently addressed the restoration 
of speckled images using TV regularization include yj, 1121 . 1161 , 



1171 . The commonalities and differences between our approach and 
the ones followed in those papers will be discussed after the detailed 
description of our method, since this discussion requires notation 
and concepts which will be introduced in the next section. 

2. PROBLEM FORMULATION 

Let y £ R" denote an n-pixels observed image, assumed to be a 
sample of a random image Y, the mean of which is the underlying 
reflectivity image X G R", ;.e.,E[Y] — x. Adopting a conditionally 
independent multiplicative noise model, we have 



Yi = XiNi, for i = 1, 



(1) 



where N £ R" is an image of independent and identically dis- 
tributed (iid) noise random variables with unit mean, E(A^i) = 1, 
following a common density p m ■ For A/-look fully developed speckle 
noise, pjv is a Gamma density with E[N] = 1, and aff — 1/M, i.e.. 



PN[n) 



r(M) 



(2) 



An additive noise model is obtained by taking logarithms of l[T). 
For some pixel of the image, the observation model becomes 



log Y — log X + log TV . 

G z W 

The density of the random variable W = log N is 



(3) 
(4) 

PG\z{g\z) = pw{g - z). (5) 

Under the regularization and Bayesian frameworks, the original 
image is inferred by solving a minimization problem with the form 



pw(-w) = pN{n = e )e 



T(M) 



Mm — e™Af 

e e , 



thus 



z G argminL(z), 

z 

where L(z) is the penalized minus log-likelihood, 

i(z) = -logPG|z(g|z) + A(^(z). 



= ME 



2. + 6" 



') +X(t){z) + A, 



(6) 



(7) 



(8) 



with A an irrelevant additive constant, the penalty/regularizer (neg- 
ative of the log-prior, from a the Bayesian perspective), and A the 
regularization parameter. 

In this work, we adopt the TV regularizer, that is. 



(z) = TV(z) = ^ V(AJz)2 + (ASz)2 



(9) 



where (A^ z and A^z) denote the horizontal and vertical first order 
differences at pixel s £ {1, . . . , n}, respectively. 

Each term (zs + e^'~^'') of l[8]l, corresponding to the negative 
log-likelihood, is strictly convex and coercive, thus so is their sum. 
Since the TV regularizer is also convex (though not strictly so), the 
objective function L possesses a unique minimizer |6). In terms 
of optimization, these are desirable properties that would not hold if 
we had formulated the inference in the original variables x, since the 
resulting negative log-likelihood is not convex; this was the approach 
followed in (JJ and [161 . 



3. BREGMAN/AUGMENTED LAGRANGIAN APPROACH 

There are several efficient algorithms to compute the TV regularized 
solution for a quadratic data term; for recent work, see f3l, fSl, fSl, 
1101 . 1191 , (21 1 > and other references therein. When the data term 
is not quadratic, as in the problem is more difficult and far less 
studied. Herein, we follow the split Bregman approach [lOJ which 
is composed of the following two steps: (splitting) a constrained 
problem equivalent to the original unconstrained one is formulated; 
(Bregman) this constrained problem is solved using the Bregman it- 
erative approach |20|. Before describing these two steps in detail, we 
briefly review the Bregman iterative approach. The reader is referred 
to 1101 . 1201 . for more details. 

3.1. Bregman Iterations 

Consider a constrained optimization problem of the form 

min E{:k) 

s.t. H{x) = 0, (10) 

with E and H convex, H differentiable, and minx H(pi) = . The 
so-called Bregman divergence associated with the convex function 
E is defined as 

Z)P(x,y) = i5(x)-i5(y)-(p,x-y), (11) 

where p belongs to the subgradient of E at y, i.e., 

p G 9£(y) = {u ; ^^(x) > £'(y) + (u,x -y), Vx G domS}. 

The Bregman iteration is given by 

x'''+^ = argminDP'°(x,x''') +//(x) (12) 

X 

= argmin£;(x) - (p'=,x-x'=) + /i'(x), (13) 

X 

where p*" G dE{x''). It has been shown that this procedure con- 
verges to a solution of l[TOj 1 10|,|20|. 

Concerning the update of p'', we have from l ll2b . that G 
9(-D^(x,x'') + fiH{x)), when this sub-differential is evaluated at 
x''+\ that is 

G a(DP(x'=+\x'=) +i?(x'=+')). 

Since it was assumed that H is differentiable, and since p'^^^ G 
dE{x'') at this point, p''^^ should be chosen as 



k + l 



= p*-' - V_H'(x 



(14) 



In the particular case where i/(x) = (r/2) || Ax - b||i, it can 
be shown (see ||20J ) that the iteration l |13t is equivalent to 

x''+^ = argmin£;(x) + -||Ax - b'-'lla (15) 
X 2 

(16) 



k+i 



b + b'= - Ax\ 



3.2. Splitting the Problem into a Constrained Formulation 

The original unconstrained problem ^ is equivalent to the con- 
strained formulation 



(z, u) = arg min L(z, u) 



s.t. 



I|z-u||^ =0, 



(17) 
(18) 



with 

n 

L(z,u) + e«=-"=) +ATV(u). (19) 

s = l 

Notice how the original variable (image) z is split into a pair of 
variables (z, u), which are decoupled in the objective function l |19l l. 

3.3. Applying Bregman Iterations 

Notice that the problem il7t-lll8t has exactly the form l IlQt . with 

X = [z^u^]^, £'(x) = L(z, u), and H{x) = (r/2)|[Ax - b||i, 
with A — [I, —I] and b = 0. Using this equivalence, the Bregman 
iteration I ll5t -(ll6t becomes 

(z''+\u''+^) = argminL(z, u) + ^llz - u - b*" 11^,(20) 

Z,U 2 

bfc+i ^ b*-(z'=-u'^). (21) 

We address the minimization in i20\ using an alternating min- 
imization scheme with respect to u and z. The complete resulting 
algorithm is summarized in Algorithm 1. 



Algorithm 1 TV restoration of multilook images. 
Initialization: z = 0, u = 0, b = 0, A, r, fc := 1. 
1: repeat 

2: for t = 1 : tm. do 

3: z^- argmin. ^^=1 + e«=-^=) + 5^||z-u'=-b'=f 

4: u'' — argminu i||u- z'' +b'=|p + 7 TV(u). 
5: end for 

6: b'=+' — b'= - (z'^ - u'=) 
7: k ■- k + 1 

8: until ||z'= - z'=-i ||i/||z'=-^ \\l < 10-" 



The minimization with respect to z, in line 3, has closed form 
in terms of the Lambert W function [7J. However, we found that the 
Newton method yields a faster solution by running just four itera- 
tions. Notice that the minimization in line 3 is in fact a set of n de- 
coupled scalar minimizations. For the minimization with respect to 
u (line 4), which is a TV denoising problem, we run a few iterations 
(typically 10) of ChamboUe's algorithm |3|. The number of inner 
iterations tm was set to one in all the experiments reported below. 
The stopping criterion (line 8) is the same as in U2J . The estimate of 
X produced by the algorithm is naturally x = e , component-wise. 

Notice how the split Bregman approach converted a difficult 
problem involving a non-quadratic term and a TV regularizer into 
two simpler problems: a decoupled minimization problem (line 3) 
and a TV denoising problem with a quadratic data term (line 4). 

3.4. Remarks 

In the case of linear constraints, the Bregman iterative procedure de- 
fined in l |13l l is equivalent to an augmented Lagrangian method 1 13 1; 
see 1181 , 1201 for proofs. It is known that the augmented Lagrangian 
is better conditioned that the standard Lagrangian for the same prob- 
lem, thus a better numerical behavior is expectable. 

TV-based image restoration under multiplicative noise was re- 
cently addressed in 1 17|. The authors apply an inverse scale space 
flow, which converges to the solution of the constrained problem of 
minimizing TV(z) under an equality constraint on the log-likelihood; 



Table 1. Experimental results (Iter denotes the number of iterations; 
Cam. is the Cameraman image). 
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this requires a carefully chosen stopping criterion, because the solu- 
tion of this constrained problem is not a good estimate. 

In 1121 . a splitting of the variable is also used to obtain an objec- 
tive function with the form 

£(z,u) = L(z,u)+a||z-u||i; (22) 

this is the so-called splitting-and-penalty method. Notice that the 
minimizers of E(z, u) converge to those of (I17l>-l ll8b only when a 
approaches infinity. However, since E{z, u) becomes severely ill- 
conditioned when Q is very large, causing numerical difficulties, it 
is only practical to minimize E{z, u) with moderate values of a; 
consequently, the solutions obtained are not minima of the regular- 
ized negative log-likelihood ([8}. 

4. EXPERIMENTS 

In this section we report experimental results comparing the perfor- 
mance of the proposed approach with that of the recent state-of-the- 
art methods in (T| and 1121 . All the experiments use synthetic data, 
in the sense that the observed image is generated according to (TJ- 
([2}, where x is a clean image. As in 1 12 1, we select the regularization 
parameter A by searching for the value leading to the lowest mean 
squared error with respect to the true image. The algorithm is ini- 
tialized with the observed noisy image. The quality of the estimates 
is assessed using the relative error (as in |T2|), 



Table[T]reports the results obtained using Lena and the Camera- 
man as original images, for the same values of the number of looks 
(A/ in ^) as used in 1 12 1. In these experiments, our method always 
achieves lower relative errors with fewer iterations, when compared 
with the methods from 1 12| and |T | (the results concerning the algo- 
rithm from 1 1 1 are those reported in 1121 ). It's important to point out 
that the computational cost of each iteration of the algorithm of |T2 | 
is essentially the same as that of our algorithm. 

Figure [T] shows the noisy and restored images, for the same ex- 
periments reported in Table [T] Finally, Figure |2] plots the evolu- 
tion of the objective function i(z*) and of the constraint function 
ll^fe _ u'°||2 along the iterations, for the example with the Cam- 
eraman image and M = 3. Observe the extremely low value of 
Hz*" — u'' II2 at the final iterations, showing that, for all practical pur- 
poses, the constraint (TS} is satisfied. 

5. CONCLUDING REMARKS 

We have proposed an approach to total variation denoising of images 
contaminated by multiplicative noise, by exploiting a split Bregman 




Fig. 1. Left column: observed noisy images. Right column: image 
estimates. First and second rows: Lena, M = 5 and M = 33. Third 
and fourth rows: Cameraman, M = 3 and M = 13. 
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Fig. 2. Evolution of the objective function L{7.^) and of the con- 
straint function Hz*^ — u''|[2, along the iterations of the algorithm, 
for the experiment with the Cameraman image and M — 3. 



technique. The proposed algorithm is very simple and, in the ex- 
periments herein reported, exhibited state of the art performance and 
speed. We are currently working on extending our methods to prob- 
lems involving linear observation operators (e.g., blur) and other re- 
lated noise models, such as Poisson. 
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